home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / minpack / dogleg.f next >
Text File  |  1996-07-19  |  5KB  |  178 lines
  1.       SUBROUTINE DOGLEG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2)
  2.       INTEGER N,LR
  3.       DOUBLE PRECISION DELTA
  4.       DOUBLE PRECISION R(LR),DIAG(N),QTB(N),X(N),WA1(N),WA2(N)
  5. C     **********
  6. C
  7. C     SUBROUTINE DOGLEG
  8. C
  9. C     GIVEN AN M BY N MATRIX A, AN N BY N NONSINGULAR DIAGONAL
  10. C     MATRIX D, AN M-VECTOR B, AND A POSITIVE NUMBER DELTA, THE
  11. C     PROBLEM IS TO DETERMINE THE CONVEX COMBINATION X OF THE
  12. C     GAUSS-NEWTON AND SCALED GRADIENT DIRECTIONS THAT MINIMIZES
  13. C     (A*X - B) IN THE LEAST SQUARES SENSE, SUBJECT TO THE
  14. C     RESTRICTION THAT THE EUCLIDEAN NORM OF D*X BE AT MOST DELTA.
  15. C
  16. C     THIS SUBROUTINE COMPLETES THE SOLUTION OF THE PROBLEM
  17. C     IF IT IS PROVIDED WITH THE NECESSARY INFORMATION FROM THE
  18. C     QR FACTORIZATION OF A. THAT IS, IF A = Q*R, WHERE Q HAS
  19. C     ORTHOGONAL COLUMNS AND R IS AN UPPER TRIANGULAR MATRIX,
  20. C     THEN DOGLEG EXPECTS THE FULL UPPER TRIANGLE OF R AND
  21. C     THE FIRST N COMPONENTS OF (Q TRANSPOSE)*B.
  22. C
  23. C     THE SUBROUTINE STATEMENT IS
  24. C
  25. C       SUBROUTINE DOGLEG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2)
  26. C
  27. C     WHERE
  28. C
  29. C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF R.
  30. C
  31. C       R IS AN INPUT ARRAY OF LENGTH LR WHICH MUST CONTAIN THE UPPER
  32. C         TRIANGULAR MATRIX R STORED BY ROWS.
  33. C
  34. C       LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
  35. C         (N*(N+1))/2.
  36. C
  37. C       DIAG IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE
  38. C         DIAGONAL ELEMENTS OF THE MATRIX D.
  39. C
  40. C       QTB IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE FIRST
  41. C         N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*B.
  42. C
  43. C       DELTA IS A POSITIVE INPUT VARIABLE WHICH SPECIFIES AN UPPER
  44. C         BOUND ON THE EUCLIDEAN NORM OF D*X.
  45. C
  46. C       X IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE DESIRED
  47. C         CONVEX COMBINATION OF THE GAUSS-NEWTON DIRECTION AND THE
  48. C         SCALED GRADIENT DIRECTION.
  49. C
  50. C       WA1 AND WA2 ARE WORK ARRAYS OF LENGTH N.
  51. C
  52. C     SUBPROGRAMS CALLED
  53. C
  54. C       MINPACK-SUPPLIED ... DPMPAR,ENORM
  55. C
  56. C       FORTRAN-SUPPLIED ... DABS,DMAX1,DMIN1,DSQRT
  57. C
  58. C     MINPACK. VERSION OF JULY 1978.
  59. C     BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
  60. C
  61. C     **********
  62.       INTEGER I,J,JJ,JP1,K,L
  63.       DOUBLE PRECISION ALPHA,BNORM,EPSMCH,GNORM,ONE,QNORM,SGNORM,SUM,
  64.      *                 TEMP,ZERO
  65.       DOUBLE PRECISION DPMPAR,ENORM
  66.       DATA ONE,ZERO /1.0D0,0.0D0/
  67. C
  68. C     EPSMCH IS THE MACHINE PRECISION.
  69. C
  70.       EPSMCH = DPMPAR(1)
  71. C
  72. C     FIRST, CALCULATE THE GAUSS-NEWTON DIRECTION.
  73. C
  74.       JJ = (N*(N + 1))/2 + 1
  75.       DO 50 K = 1, N
  76.          J = N - K + 1
  77.          JP1 = J + 1
  78.          JJ = JJ - K
  79.          L = JJ + 1
  80.          SUM = ZERO
  81.          IF (N .LT. JP1) GO TO 20
  82.          DO 10 I = JP1, N
  83.             SUM = SUM + R(L)*X(I)
  84.             L = L + 1
  85.    10       CONTINUE
  86.    20    CONTINUE
  87.          TEMP = R(JJ)
  88.          IF (TEMP .NE. ZERO) GO TO 40
  89.          L = J
  90.          DO 30 I = 1, J
  91.             TEMP = DMAX1(TEMP,DABS(R(L)))
  92.             L = L + N - I
  93.    30       CONTINUE
  94.          TEMP = EPSMCH*TEMP
  95.          IF (TEMP .EQ. ZERO) TEMP = EPSMCH
  96.    40    CONTINUE
  97.          X(J) = (QTB(J) - SUM)/TEMP
  98.    50    CONTINUE
  99. C
  100. C     TEST WHETHER THE GAUSS-NEWTON DIRECTION IS ACCEPTABLE.
  101. C
  102.       DO 60 J = 1, N
  103.          WA1(J) = ZERO
  104.          WA2(J) = DIAG(J)*X(J)
  105.    60    CONTINUE
  106.       QNORM = ENORM(N,WA2)
  107.       IF (QNORM .LE. DELTA) GO TO 140
  108. C
  109. C     THE GAUSS-NEWTON DIRECTION IS NOT ACCEPTABLE.
  110. C     NEXT, CALCULATE THE SCALED GRADIENT DIRECTION.
  111. C
  112.       L = 1
  113.       DO 80 J = 1, N
  114.          TEMP = QTB(J)
  115.          DO 70 I = J, N
  116.             WA1(I) = WA1(I) + R(L)*TEMP
  117.             L = L + 1
  118.    70       CONTINUE
  119.          WA1(J) = WA1(J)/DIAG(J)
  120.    80    CONTINUE
  121. C
  122. C     CALCULATE THE NORM OF THE SCALED GRADIENT DIRECTION,
  123. C     NORMALIZE, AND RESCALE THE GRADIENT.
  124. C
  125.       GNORM = ENORM(N,WA1)
  126.       SGNORM = ZERO
  127.       ALPHA = DELTA/QNORM
  128.       IF (GNORM .EQ. ZERO) GO TO 120
  129.       DO 90 J = 1, N
  130.          WA1(J) = (WA1(J)/GNORM)/DIAG(J)
  131.    90    CONTINUE
  132. C
  133. C     CALCULATE THE POINT ALONG THE SCALED GRADIENT
  134. C     AT WHICH THE QUADRATIC IS MINIMIZED.
  135. C
  136.       L = 1
  137.       DO 110 J = 1, N
  138.          SUM = ZERO
  139.          DO 100 I = J, N
  140.             SUM = SUM + R(L)*WA1(I)
  141.             L = L + 1
  142.   100       CONTINUE
  143.          WA2(J) = SUM
  144.   110    CONTINUE
  145.       TEMP = ENORM(N,WA2)
  146.       SGNORM = (GNORM/TEMP)/TEMP
  147. C
  148. C     TEST WHETHER THE SCALED GRADIENT DIRECTION IS ACCEPTABLE.
  149. C
  150.       ALPHA = ZERO
  151.       IF (SGNORM .GE. DELTA) GO TO 120
  152. C
  153. C     THE SCALED GRADIENT DIRECTION IS NOT ACCEPTABLE.
  154. C     FINALLY, CALCULATE THE POINT ALONG THE DOGLEG
  155. C     AT WHICH THE QUADRATIC IS MINIMIZED.
  156. C
  157.       BNORM = ENORM(N,QTB)
  158.       TEMP = (BNORM/GNORM)*(BNORM/QNORM)*(SGNORM/DELTA)
  159.       TEMP = TEMP - (DELTA/QNORM)*(SGNORM/DELTA)**2
  160.      *       + DSQRT((TEMP-(DELTA/QNORM))**2
  161.      *               +(ONE-(DELTA/QNORM)**2)*(ONE-(SGNORM/DELTA)**2))
  162.       ALPHA = ((DELTA/QNORM)*(ONE - (SGNORM/DELTA)**2))/TEMP
  163.   120 CONTINUE
  164. C
  165. C     FORM APPROPRIATE CONVEX COMBINATION OF THE GAUSS-NEWTON
  166. C     DIRECTION AND THE SCALED GRADIENT DIRECTION.
  167. C
  168.       TEMP = (ONE - ALPHA)*DMIN1(SGNORM,DELTA)
  169.       DO 130 J = 1, N
  170.          X(J) = TEMP*WA1(J) + ALPHA*X(J)
  171.   130    CONTINUE
  172.   140 CONTINUE
  173.       RETURN
  174. C
  175. C     LAST CARD OF SUBROUTINE DOGLEG.
  176. C
  177.       END
  178.